######################################################################################
# REPLICATION CODE:
# AUTHORS: Sebastian Juhl and Laron K. Williams
# ARTICLE: "Learning at Home and Abroad: How Competition Conditions the
#           Diffusion of Party Strategies"
# JOURNAL: British Journal of Political Science
######################################################################################

# NOTE:
# This code produces all figures and tables reported in the main
# article and the supplementary materials. The results will be stored
# in the subfolder "Output" which will be created unless it exists already.

# set working directory
setwd(YOUR_DIRECTORY)

# clear working environment
rm(list=ls(all=T))

# create empty output folder
ifelse(!dir.exists(file.path(paste0(getwd(),"/Output")))
       ,dir.create(file.path(paste0(getwd(),"/Output"))), FALSE)

### PACKAGES
library(MASS)
library(dplyr)
library(ggplot2)

### LOAD DATASETS
# load data
load("./Data/Juhl_Williams_Data.RData")
# load connectivity matrices
load("./Data/Matrices.RData")

# select complete cases
complete <- complete.cases(data)
cmp <- cmp[complete,]
data <- data[complete,]


# define spatial lags
W.f.win <- W.f.win[complete,complete]
W.f.lose <- W.f.lose[complete,complete]
W.d.win <- W.d.win[complete,complete]
W.d.lose <- W.d.lose[complete,complete]
data$f.bloc.win <- W.f.win %*% data$issue
data$f.bloc.lose <- W.f.lose %*% data$issue
data$d.bloc.win.lag <- W.d.win %*% data$issue
data$d.bloc.lose.lag <- W.d.lose %*% data$issue


# sum of links
s.f.bloc.win <- apply(W.f.win,1,sum)
s.f.bloc.lose <- apply(W.f.lose,1,sum)
s.d.bloc.win.lag <- apply(W.d.win,1,sum)
s.d.bloc.lose.lag <- apply(W.d.lose,1,sum)

# average number of neighbors
avg.f.bloc.win <- mean(apply(W.f.win,1,sum))
avg.f.bloc.lose <- mean(apply(W.f.lose,1,sum))
avg.d.bloc.win <- mean(apply(W.d.win,1,sum))
avg.d.bloc.lose <- mean(apply(W.d.lose,1,sum))


###################
# Table 2
###################
nams <- c("y","DV_t-1","W^Fwy","W^Fly","W^Dwy","W^Dly"
          ,"VSGreen_t-1","AvgGreen_n-1","GDP_n-1","Incumbent_t"
          ,"Incumbent_DVt-1","VS_t-1")
tabs <- ifelse(nchar(nams)<=3,"\t\t\t\t"
               ,ifelse(nchar(nams)<=6,"\t\t\t"
               ,ifelse(nchar(nams)<=9,"\t\t"
               ,ifelse(nchar(nams)<=12,"\t",""))))
pos <- c(1:2,10:13,3,9,5:6,8,7)
Mean <- sprintf("%.3f",round(apply(data[,pos],2,mean),3))
SD <- sprintf("%.3f",round(apply(data[,pos],2,sd),3))
Min <- round(apply(data[,pos],2,min),3)
Max <- round(apply(data[,pos],2,max),3)

cat("TABLE 2: Descriptive Statistics\n\n"
    ,append=F, file="./Output/Tab2.txt")
cat("____________________________________________________\n"
    ,append=T, file="./Output/Tab2.txt")
for(i in 1:length(nams)){
  cat(nams[i],tabs[i],Mean[i],"\t\t",SD[i],"\t\t",Min[i],"\t\t",Max[i],"\n"
      ,append=T, file="./Output/Tab2.txt")
}


###################
# Table 3
###################
# Model 1
m1 <- glm(issue ~ issue.lag
          + f.bloc.win + f.bloc.lose
          + vsgreen.lag
          + vsgreen.lag2
          + globalgreen.lag
          + gdp.lag
          + incumbent
          + incumbent.emphasis.lag
          + prevvote
          + factor(cmp$party)
          + factor(cmp$year)
          ,data=data
          ,family=quasibinomial("logit"))
c1 <- sprintf("%.3f",round(coef(m1),3)[1:11])
s1 <- sprintf("%.3f",round(summary(m1)$coefficients[,"Std. Error"],3)[1:11])
pos1 <- c(1:4,7:13)
coef1 <- rep("-",length(c1))
coef1[pos1] <- c1
se1 <- rep("",length(s1))
se1[pos1] <- s1

# Model 2
m2 <- glm(issue ~ issue.lag
          + d.bloc.win.lag + d.bloc.lose.lag
          + vsgreen.lag
          + vsgreen.lag2
          + globalgreen.lag
          + gdp.lag
          + incumbent
          + incumbent.emphasis.lag
          + prevvote
          + factor(cmp$party)
          + factor(cmp$year)
          ,data=data
          ,family=quasibinomial("logit"))
c2 <- sprintf("%.3f",round(coef(m2),3)[1:11])
s2 <- sprintf("%.3f",round(summary(m2)$coefficients[,"Std. Error"],3)[1:11])
pos2 <- c(1:2,5:13)
coef2 <- rep("-",length(c2))
coef2[pos2] <- c2
se2 <- rep("",length(s2))
se2[pos2] <- s2

# Model 3
m3 <- glm(issue ~ issue.lag
          + f.bloc.win + f.bloc.lose
          + d.bloc.win.lag + d.bloc.lose.lag
          + vsgreen.lag
          + vsgreen.lag2
          + globalgreen.lag
          + gdp.lag
          + incumbent
          + incumbent.emphasis.lag
          + prevvote
          + factor(cmp$party)
          + factor(cmp$year)
          ,data=data
          ,family=quasibinomial("logit"))
coef3 <- sprintf("%.3f",round(coef(m3),3)[1:13])
se3 <- sprintf("%.3f",round(summary(m3)$coefficients[,"Std. Error"],3)[1:13])


# phi
phi <- c(round(summary(m1)$dispersion,4)
         ,round(summary(m2)$dispersion,4)
         ,round(summary(m3)$dispersion,4))
# Observations
obs <- c(length(m1$residuals),length(m2$residuals),length(m3$residuals))
# RMSE
RMSE <- c(round(sqrt(mean((m1$fitted.values - m1$y)^2)),4)
          ,round(sqrt(mean((m2$fitted.values - m2$y)^2)),4)
          ,round(sqrt(mean((m3$fitted.values - m3$y)^2)),4))

cat("TABLE 3: Fractional Logit Model Estimates of Environmental Issues Emphasis\n\n"
    ,append=F, file="./Output/Tab3.txt")
cat("\t\t\t\t Model 1 \t\t Model 2 \t\t Model 3\n"
    ,append=T, file="./Output/Tab3.txt")
cat("__________________________________________________________\n"
    ,append=T, file="./Output/Tab3.txt")
cat("Constant\t\t",coef1[1],"\t\t",coef2[1],"\t\t",coef3[1],"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("\t\t\t\t",paste0("(",se1[1],")"),"\t\t",paste0("(",se2[1],")"),"\t\t",paste0("(",se3[1],")"),"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("DV_t-1\t\t\t",coef1[2],"\t\t",coef2[2],"\t\t",coef3[2],"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("\t\t\t\t",paste0("(",se1[2],")"),"\t\t",paste0("(",se2[2],")"),"\t\t",paste0("(",se3[2],")"),"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("theta_1W^Fwy\t",coef1[3],"\t\t\t",coef2[3],"\t\t\t\t",coef3[3],"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("\t\t\t\t",paste0("(",se1[3],")"),"\t\t",paste0("(",se2[3],")"),"\t\t\t",paste0("(",se3[3],")"),"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("theta_2W^Fly\t",coef1[4],"\t\t",coef2[4],"\t\t\t\t",coef3[4],"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("\t\t\t\t",paste0("(",se1[4],")"),"\t\t",paste0("(",se2[4],")"),"\t\t\t",paste0("(",se3[4],")"),"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("theta_3W^Dwy\t",coef1[5],"\t\t\t\t",coef2[5],"\t\t",coef3[5],"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("\t\t\t\t",paste0("(",se1[5],")"),"\t\t\t",paste0("(",se2[5],")"),"\t\t",paste0("(",se3[5],")"),"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("theta_4W^Dly\t",coef1[6],"\t\t\t\t",coef2[6],"\t\t",coef3[6],"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("\t\t\t\t",paste0("(",se1[6],")"),"\t\t\t",paste0("(",se2[6],")"),"\t\t",paste0("(",se3[6],")"),"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("VSGreen_t-1\t\t",coef1[7],"\t\t",coef2[7],"\t\t",coef3[7],"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("\t\t\t\t",paste0("(",se1[7],")"),"\t\t",paste0("(",se2[7],")"),"\t\t",paste0("(",se3[7],")"),"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("VSGreen^2_t-1\t",coef1[8],"\t\t\t",coef2[8],"\t\t\t",coef3[8],"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("\t\t\t\t",paste0("(",se1[8],")"),"\t\t",paste0("(",se2[8],")"),"\t\t",paste0("(",se3[8],")"),"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("AvgGreen_n-1\t",coef1[9],"\t\t\t",coef2[9],"\t\t\t",coef3[9],"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("\t\t\t\t",paste0("(",se1[9],")"),"\t\t",paste0("(",se2[9],")"),"\t\t",paste0("(",se3[9],")"),"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("GDP_n-1\t\t\t",coef1[10],"\t\t\t",coef2[10],"\t\t\t",coef3[10],"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("\t\t\t\t",paste0("(",se1[10],")"),"\t\t",paste0("(",se2[10],")"),"\t\t",paste0("(",se3[10],")"),"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("Incumbent_t\t\t",coef1[11],"\t\t",coef2[11],"\t\t",coef3[11],"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("\t\t\t\t",paste0("(",se1[11],")"),"\t\t",paste0("(",se2[11],")"),"\t\t",paste0("(",se3[11],")"),"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("Incumbent_DV_t-1",coef1[12],"\t\t\t",coef2[12],"\t\t\t",coef3[12],"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("\t\t\t\t",paste0("(",se1[12],")"),"\t\t",paste0("(",se2[12],")"),"\t\t",paste0("(",se3[12],")"),"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("VS_t-1\t\t\t",coef1[13],"\t\t\t",coef2[13],"\t\t\t",coef3[13],"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("\t\t\t\t",paste0("(",se1[13],")"),"\t\t",paste0("(",se2[13],")"),"\t\t",paste0("(",se3[13],")"),"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("Party and\t\t","Yes \t\t\t Yes \t\t\t Yes\n"
    ,append=T, file="./Output/Tab3.txt")
cat("Year FEs\n"
    ,append=T, file="./Output/Tab3.txt")
cat("__________________________________________________________\n"
    ,append=T, file="./Output/Tab3.txt")
cat("phi\t\t\t\t",phi[1],"\t\t",phi[2],"\t\t",phi[3],"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("Observations\t",obs[1],"\t\t\t",obs[2],"\t\t\t",obs[3],"\n"
    ,append=T, file="./Output/Tab3.txt")
cat("RMSE\t\t\t",RMSE[1],"\t\t",RMSE[2],"\t\t",RMSE[3]
    ,append=T, file="./Output/Tab3.txt")



###################
# Table C.1
###################
m5 <- glm(issue ~ issue.lag
          + f.bloc.win + f.bloc.lose
          + d.bloc.win.lag + d.bloc.lose.lag
          + I(vsgreen.lag*f.bloc.win) + I(vsgreen.lag*f.bloc.lose)
          + vsgreen.lag
          + vsgreen.lag2
          + globalgreen.lag
          + gdp.lag
          + incumbent
          + incumbent.emphasis.lag
          + prevvote
          + factor(cmp$party)
          + factor(cmp$year)
          ,data=data
          ,family=quasibinomial("logit"))
coef5 <- sprintf("%.3f",round(coef(m5),3)[1:15])
se5 <- sprintf("%.3f",round(summary(m5)$coefficients[,"Std. Error"],3)[1:15])

# phi
phi <- c(round(summary(m3)$dispersion,4),round(summary(m5)$dispersion,4))
# Observations
obs <- c(length(m3$residuals),length(m5$residuals))
# RMSE
RMSE <- c(round(sqrt(mean((m3$fitted.values - m3$y)^2)),4)
          ,round(sqrt(mean((m5$fitted.values - m5$y)^2)),4))


cat("TABLE C.1: Fractional Logit Model Estimates of Environmental Issues Emphasis with Conditional Spatial Effects\n\n"
    ,append=F, file="./Output/TabC.1.txt")
cat("\t\t\t\t Model 3 \t\t Model 5\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("___________________________________________\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("Constant\t\t",coef3[1],"\t\t",coef5[1],"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("\t\t\t\t",paste0("(",se3[1],")"),"\t\t",paste0("(",se5[1],")"),"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("DV_t-1\t\t\t",coef3[2],"\t\t",coef5[2],"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("\t\t\t\t",paste0("(",se3[2],")"),"\t\t",paste0("(",se5[2],")"),"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("theta_1W^Fwy\t",coef3[3],"\t\t\t",coef5[3],"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("\t\t\t\t",paste0("(",se3[3],")"),"\t\t",paste0("(",se5[3],")"),"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("theta_2W^Fly\t",coef3[4],"\t\t",coef5[4],"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("\t\t\t\t",paste0("(",se3[4],")"),"\t\t",paste0("(",se5[4],")"),"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("theta_3W^Dwy\t",coef3[5],"\t\t\t",coef5[5],"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("\t\t\t\t",paste0("(",se3[5],")"),"\t\t",paste0("(",se5[5],")"),"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("theta_4W^Dly\t",coef3[6],"\t\t",coef5[6],"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("\t\t\t\t",paste0("(",se3[6],")"),"\t\t",paste0("(",se5[6],")"),"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("VSGreen_t-1 x\t","-","\t\t\t\t",coef5[7],"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("W^Fwy\t\t\t","()","\t\t\t",paste0("(",se5[7],")"),"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("VSGreen_t-1 x\t","-","\t\t\t\t",coef5[8],"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("W^Fly\t\t\t","()","\t\t\t",paste0("(",se5[8],")"),"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("VSGreen_t-1\t\t",coef3[7],"\t\t",coef5[9],"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("\t\t\t\t",paste0("(",se3[7],")"),"\t\t",paste0("(",se5[9],")"),"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("VSGreen^2_t-1\t",coef3[8],"\t\t\t",coef5[10],"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("\t\t\t\t",paste0("(",se3[8],")"),"\t\t",paste0("(",se5[10],")"),"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("AvgGreen_n-1\t",coef3[9],"\t\t\t",coef5[11],"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("\t\t\t\t",paste0("(",se3[9],")"),"\t\t",paste0("(",se5[11],")"),"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("GDP_n-1\t\t\t",coef3[10],"\t\t\t",coef5[12],"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("\t\t\t\t",paste0("(",se3[10],")"),"\t\t",paste0("(",se5[12],")"),"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("Incumbent_t\t\t",coef3[11],"\t\t",coef5[13],"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("\t\t\t\t",paste0("(",se3[11],")"),"\t\t",paste0("(",se5[13],")"),"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("Incumbent_DV_t-1",coef3[12],"\t\t\t",coef5[14],"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("\t\t\t\t",paste0("(",se3[12],")"),"\t\t",paste0("(",se5[14],")"),"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("VS_t-1\t\t\t",coef3[13],"\t\t\t",coef5[15],"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("\t\t\t\t",paste0("(",se3[13],")"),"\t\t",paste0("(",se5[15],")"),"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("Party and\t\t","Yes \t\t\t Yes\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("Year FEs\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("___________________________________________\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("phi\t\t\t\t",phi[1],"\t\t",phi[2],"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("Observations\t",obs[1],"\t\t\t",obs[2],"\n"
    ,append=T, file="./Output/TabC.1.txt")
cat("RMSE\t\t\t",RMSE[1],"\t\t",RMSE[2]
    ,append=T, file="./Output/TabC.1.txt")



###################
# Table D.1
###################
# define spatial lags
W.d.all <- W.d.all[complete,complete]
data$d.all.lag <- W.d.all %*% data$issue

# Model 6
m6 <- glm(issue ~ issue.lag
          + f.bloc.win + f.bloc.lose
          + d.all.lag
          + vsgreen.lag
          + vsgreen.lag2
          + globalgreen.lag
          + gdp.lag
          + incumbent
          + incumbent.emphasis.lag
          + prevvote
          + factor(cmp$party)
          + factor(cmp$year)
          ,data=data
          ,family=quasibinomial("logit"))
c6 <- sprintf("%.3f",round(coef(m6),3)[1:14])
s6 <- sprintf("%.3f",round(summary(m6)$coefficients[,"Std. Error"],3)[1:14])
pos6 <- c(1:4,7:14)
coef6 <- rep("-",length(c6))
coef6[pos6] <- c6
se6 <- rep("",length(s6))
se6[pos6] <- s6

# Model 7
m7 <- glm(issue ~ issue.lag
          + f.bloc.win + f.bloc.lose
          + d.bloc.win.lag + d.bloc.lose.lag
          + d.all.lag
          + vsgreen.lag
          + vsgreen.lag2
          + globalgreen.lag
          + gdp.lag
          + incumbent
          + incumbent.emphasis.lag
          + prevvote
          + factor(cmp$party)
          + factor(cmp$year)
          ,data=data
          ,family=quasibinomial("logit"))
coef7 <- sprintf("%.3f",round(coef(m7),3)[1:14])
se7 <- sprintf("%.3f",round(summary(m7)$coefficients[,"Std. Error"],3)[1:14])

# phi
phi <- c(round(summary(m3)$dispersion,4)
         ,round(summary(m6)$dispersion,4)
         ,round(summary(m7)$dispersion,4))
# Observations
obs <- c(length(m3$residuals),length(m6$residuals),length(m7$residuals))
# RMSE
RMSE <- c(round(sqrt(mean((m3$fitted.values - m3$y)^2)),4)
          ,round(sqrt(mean((m6$fitted.values - m6$y)^2)),4)
          ,round(sqrt(mean((m7$fitted.values - m7$y)^2)),4))


cat("TABLE D.1: Fractional Logit Model Estimates of Environmental Issues Emphasis with Spatial Lags Connecting all Domestic Parties\n\n"
    ,append=F, file="./Output/TabD.1.txt")
cat("\t\t\t\t Model 3 \t\t Model 6 \t\t Model 7\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("__________________________________________________________\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("Constant\t\t",coef3[1],"\t\t",coef6[1],"\t\t",coef7[1],"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("\t\t\t\t",paste0("(",se3[1],")"),"\t\t",paste0("(",se6[1],")"),"\t\t",paste0("(",se7[1],")"),"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("DV_t-1\t\t\t",coef3[2],"\t\t",coef6[2],"\t\t",coef7[2],"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("\t\t\t\t",paste0("(",se3[2],")"),"\t\t",paste0("(",se6[2],")"),"\t\t",paste0("(",se7[2],")"),"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("theta_1W^Fwy\t",coef3[3],"\t\t\t",coef6[3],"\t\t\t",coef7[3],"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("\t\t\t\t",paste0("(",se3[3],")"),"\t\t",paste0("(",se6[3],")"),"\t\t",paste0("(",se7[3],")"),"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("theta_2W^Fly\t",coef3[4],"\t\t",coef6[4],"\t\t",coef7[4],"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("\t\t\t\t",paste0("(",se3[4],")"),"\t\t",paste0("(",se6[4],")"),"\t\t",paste0("(",se7[4],")"),"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("theta_3W^Dwy\t",coef3[5],"\t\t\t",coef6[5],"\t\t\t\t",coef7[5],"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("\t\t\t\t",paste0("(",se3[5],")"),"\t\t",paste0("(",se6[5],")"),"\t\t\t",paste0("(",se7[5],")"),"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("theta_4W^Dly\t",coef3[6],"\t\t",coef6[6],"\t\t\t\t",coef7[6],"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("\t\t\t\t",paste0("(",se3[6],")"),"\t\t",paste0("(",se6[6],")"),"\t\t\t",paste0("(",se7[6],")"),"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("theta_5W^Dally\t","-","\t\t\t\t",coef6[7],"\t\t\t",coef7[7],"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("\t\t\t\t","()","\t\t\t",paste0("(",se6[7],")"),"\t\t",paste0("(",se7[7],")"),"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("VSGreen_t-1\t\t",coef3[7],"\t\t",coef6[8],"\t\t",coef7[8],"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("\t\t\t\t",paste0("(",se3[7],")"),"\t\t",paste0("(",se6[8],")"),"\t\t",paste0("(",se7[8],")"),"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("VSGreen^2_t-1\t",coef3[8],"\t\t\t",coef6[9],"\t\t\t",coef7[9],"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("\t\t\t\t",paste0("(",se3[8],")"),"\t\t",paste0("(",se6[9],")"),"\t\t",paste0("(",se7[9],")"),"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("AvgGreen_n-1\t",coef3[9],"\t\t\t",coef6[10],"\t\t\t",coef7[10],"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("\t\t\t\t",paste0("(",se3[9],")"),"\t\t",paste0("(",se6[10],")"),"\t\t",paste0("(",se7[10],")"),"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("GDP_n-1\t\t\t",coef3[10],"\t\t\t",coef6[11],"\t\t\t",coef7[11],"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("\t\t\t\t",paste0("(",se3[10],")"),"\t\t",paste0("(",se6[11],")"),"\t\t",paste0("(",se7[11],")"),"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("Incumbent_t\t\t",coef3[11],"\t\t",coef6[12],"\t\t",coef7[12],"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("\t\t\t\t",paste0("(",se3[11],")"),"\t\t",paste0("(",se6[12],")"),"\t\t",paste0("(",se7[12],")"),"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("Incumbent_DV_t-1",coef3[12],"\t\t\t",coef6[13],"\t\t\t",coef7[13],"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("\t\t\t\t",paste0("(",se3[12],")"),"\t\t",paste0("(",se6[13],")"),"\t\t",paste0("(",se7[13],")"),"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("VS_t-1\t\t\t",coef3[13],"\t\t\t",coef6[14],"\t\t\t",coef7[14],"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("\t\t\t\t",paste0("(",se3[13],")"),"\t\t",paste0("(",se6[14],")"),"\t\t",paste0("(",se7[14],")"),"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("Party and\t\t","Yes \t\t\t Yes \t\t\t Yes\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("Year FEs\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("__________________________________________________________\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("phi\t\t\t\t",phi[1],"\t\t",phi[2],"\t\t",phi[3],"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("Observations\t",obs[1],"\t\t\t",obs[2],"\t\t\t",obs[3],"\n"
    ,append=T, file="./Output/TabD.1.txt")
cat("RMSE\t\t\t",RMSE[1],"\t\t",RMSE[2],"\t\t",RMSE[3]
    ,append=T, file="./Output/TabD.1.txt")



###################
# Table E.1
###################
vcov <- vcov(m3)
vcov <- vcov[rownames(vcov)=="vsgreen.lag" | rownames(vcov)=="vsgreen.lag2",colnames(vcov)=="vsgreen.lag" | colnames(vcov)=="vsgreen.lag2"]
coefs <- m3$coefficients[names(m3$coefficients)=="vsgreen.lag" | names(m3$coefficients)=="vsgreen.lag2"]
xl <- min(data$vsgreen.lag)
xh <- max(data$vsgreen.lag)
# slopes and SEs
# xl
slope_xl <- coefs[1] + 2*coefs[2] * xl
se_slope_xl <- sqrt((vcov[1,1] + 4*xl*vcov[1,2] + 4*xl^2*vcov[2,2]))
#xh
slope_xh <- coefs[1] + 2*coefs[2] * xh
se_slope_xh <- sqrt((vcov[1,1] + 4*xh*vcov[1,2] + 4*xh^2*vcov[2,2]))
# test statistic
tl <- (coefs[1] + 2*coefs[2]*xl)/(sqrt(vcov[1,1]+4*xl^2*vcov[2,2]+4*xl*vcov[1,2]))
th <- (coefs[1] + 2*coefs[2]*xh)/(sqrt(vcov[1,1]+4*xh^2*vcov[2,2]+4*xh*vcov[1,2]))

cat("TABLE E.1: Test for an U-shaped Relationship\n\n"
    ,append=F, file="./Output/TabE.1.txt")
cat("_____________________________________________\n"
    ,append=T, file="./Output/TabE.1.txt")
cat("\t\t\t\t X_min=0 \t\t X_max=14.35\n"
    ,append=T, file="./Output/TabE.1.txt")
cat("_____________________________________________\n"
    ,append=T, file="./Output/TabE.1.txt")
cat("Slope\t\t\t",sprintf("%.3f",slope_xl),"\t\t"
    ,sprintf("%.3f",slope_xh),"\n"
    ,append=T, file="./Output/TabE.1.txt")
cat("\t\t\t\t",paste0("(",sprintf("%.3f",se_slope_xl),")"),"\t\t"
    ,paste0("(",sprintf("%.3f",se_slope_xh),")"),"\n"
    ,append=T, file="./Output/TabE.1.txt")
cat("Lind-Mehlum\t\t",sprintf("%.3f",round(tl,3)),"\t\t"
    ,sprintf("%.3f",round(th,3)),"\n"
    ,append=T, file="./Output/TabE.1.txt")
cat("test statistic\n"
    ,append=T, file="./Output/TabE.1.txt")
cat("p-value\t\t\t",sprintf("%.3f",round(2*pt(abs(tl),(nrow(data)-1),lower=F),3)),"\t\t\t"
    ,sprintf("%.3f",round(2*pt(abs(th),(nrow(data)-1),lower=F),3)),"\n"
    ,append=T, file="./Output/TabE.1.txt")
cat("_____________________________________________"
    ,append=T, file="./Output/TabE.1.txt")



### Average Marginal Effects
set.seed(123)
# Extract information for Simulations
coefs <- m3$coefficients[1:13]
vcov <- vcov(m3)[1:13,1:13]
nsim <- 10000
sims <- mvrnorm(nsim,mu=coefs,Sigma=vcov)
ame.sim <- matrix(NA,nrow=nsim,ncol=4)
colnames(ame.sim) <- c("f.bloc.win","f.bloc.lose"
                       ,"d.bloc.win.lag","d.bloc.lose.lag")
greenvote <- seq(from=0, to=17,by=.1)
gvote <- matrix(NA,nrow=nsim,ncol=length(greenvote))
for(i in 1:nsim){
  # AME of spatial coefficients
  scenario <- cbind(1,data$issue.lag
                    ,data$f.bloc.win
                    ,data$f.bloc.lose
                    ,data$d.bloc.win.lag
                    ,data$d.bloc.lose.lag
                    ,data$vsgreen.lag,data$vsgreen.lag2
                    ,data$globalgreen.lag
                    ,data$gdp.lag,data$incumbent
                    ,data$incumbent.emphasis.lag
                    ,data$prevvote)
  # xb
  xb <- scenario %*% sims[i,]
  # get probabilities
  p <- exp(xb)/(1+exp(xb))
  
  # AMEs spatial short-term effects
  ame.sim[i,"f.bloc.win"] <- 1/nrow(data) * sum((sims[i,"f.bloc.win"]*s.f.bloc.win)*p*(1-p))
  ame.sim[i,"f.bloc.lose"] <- 1/nrow(data) * sum((sims[i,"f.bloc.lose"]*s.f.bloc.lose)*p*(1-p))
  ame.sim[i,"d.bloc.win.lag"] <- 1/nrow(data) * sum((sims[i,"d.bloc.win.lag"]*s.d.bloc.win.lag)*p*(1-p))
  ame.sim[i,"d.bloc.lose.lag"] <- 1/nrow(data) * sum((sims[i,"d.bloc.lose.lag"]*s.d.bloc.lose.lag)*p*(1-p))
  
  # Marginal Effect of Green Vote Share
  for(g in 1:ncol(gvote)){
    gvote[i,g] <- 1/nrow(data) * sum((sims[i,"vsgreen.lag"] + 2*sims[i,"vsgreen.lag2"] * greenvote[g]) *p*(1-p))
  }
}

ame <- matrix(NA,nrow=4,ncol=3)
colnames(ame) <- c("lowerCI","Estimate","upperCI")
rownames(ame) <- c("f.bloc.win","f.bloc.lose","d.bloc.win.lag","d.bloc.lose.lag")
ame[,"Estimate"] <- apply(ame.sim,2,mean)
ame[,"lowerCI"] <- apply(ame.sim,2,quantile, probs=.025)
ame[,"upperCI"] <- apply(ame.sim,2,quantile, probs=.975)


###################
# Figure 1
###################
# Average Marginal Effects (AME) Plot
pdf("./Output/Figure_1.pdf",width=5,height=3)
par(mar=c(3.5,5,0.1,2))
plot(x=ame[,"Estimate"],y=c(4,3,2,1),ylim=c(0,5),xlim=c(-.45,.85),pch=4,axes=F,ann=F)
lines(y=c(4,4),x=c(ame["f.bloc.win","lowerCI"],ame["f.bloc.win","upperCI"]))
lines(y=c(3,3),x=c(ame["f.bloc.lose","lowerCI"],ame["f.bloc.lose","upperCI"]))
lines(y=c(2,2),x=c(ame["d.bloc.win.lag","lowerCI"],ame["d.bloc.win.lag","upperCI"]))
lines(y=c(1,1),x=c(ame["d.bloc.lose.lag","lowerCI"],ame["d.bloc.lose.lag","upperCI"]))
text(x=c(ame[1,2],ame[2,2]-.005,ame[3,2]+.07,ame[4,2]-.05),y=c(4,3,2,1),labels=round(ame[,2],3),pos=3,cex=.7)
abline(v=0,lty=2)
axis(1)
axis(2,c(-1,1,2,3,4,6)
     ,labels=c("",expression(paste(theta[4], bold(W)^{D[L]}, bold(y)))
               ,expression(paste(theta[3], bold(W)^{D[W]}, bold(y)))
               ,expression(paste(theta[2], bold(W)^{F[L]}, bold(y)))
               ,expression(paste(theta[1], bold(W)^{F[W]}, bold(y))),"")
     ,las=1)
mtext("Effect Size",1,line=2.4)
box("plot",lwd=2)
dev.off()


###################
# Figure E.1
###################
# Marginal Effects Plot Greenvote
est <- apply(gvote,2,mean)
low <- apply(gvote,2,quantile,probs=.025)
up <- apply(gvote,2,quantile,probs=.975)
pdf("./Output/Figure_E.1.pdf",width=5,height=3)
par(mar=c(3.5,4.9,0.1,2))
plot(x=greenvote,y=est,ylim=c(-.035,.035),xlim=c(0.5,15),type="l",lwd=2.5,axes=F,ann=F)
lines(greenvote,low)
lines(greenvote,up)
abline(h=0,lty=2)
rug(data$vsgreen.lag,col="gray70")
axis(1,c(-10,0,5,10,15,20),labels=c("","0","5","10","15",""),cex.axis=.9)
axis(2,las=1,cex.axis=.9)
mtext("Green Vote Share at t-1",1,line=2.4)
mtext("Marginal Effect",2,line=3)
box("plot",lwd=2)
dev.off()


###################
# Simulation Study
###################
# define selector
sel <- cmp$countryname=="Sweden"&cmp$year==2010&cmp$parfam==60

set.seed(123)
# Extract information for Simulations
coefs <- m3$coefficients
vcov <- vcov(m3)
nsim <- 10000
sims <- mvrnorm(nsim,mu=coefs[!is.na(coefs)]
                ,Sigma=vcov[!is.na(coefs),!is.na(coefs)])

fes <- rep(0,sum(grepl("^factor\\(cmp\\$",names(coefs))))
fes[names(coefs[grepl("^factor\\(cmp\\$",names(coefs))])==paste0("factor(cmp$party)",cmp$party[sel])] <- 1
fes[names(coefs[grepl("^factor\\(cmp\\$",names(coefs))])==paste0("factor(cmp$year)",cmp$year[sel])] <- 1

# define true scenario (conservative party in Sweden 2010)
truesweden2010 <- c(1,data$issue.lag[sel]
                           ,data$f.bloc.win[sel]
                           ,data$f.bloc.lose[sel]
                           ,data$d.bloc.win.lag[sel]
                           ,data$d.bloc.lose.lag[sel]
                           ,data$vsgreen.lag[sel]
                           ,data$vsgreen.lag[sel]^2
                           ,data$globalgreen.lag[sel]
                           ,data$gdp.lag[sel]
                           ,data$incumbent[sel] # gov party
                           ,data$incumbent.emphasis.lag[sel]
                           ,data$prevvote[sel]
                           ,fes)
truesweden2010 <- truesweden2010[!is.na(coefs)]

# define fabricated y (issue emphasis of conservative party 2010 in Sweden if vsgreen.lag doubles)
fakesweden2010 <- c(1,data$issue.lag[sel]
                        ,data$f.bloc.win[sel]
                        ,data$f.bloc.lose[sel]
                        ,data$d.bloc.win.lag[sel]
                        ,data$d.bloc.lose.lag[sel]
                        ,data$vsgreen.lag[sel] *2
                        ,(data$vsgreen.lag[sel] *2)^2
                        ,data$globalgreen.lag[sel]
                        ,data$gdp.lag[sel]
                        ,data$incumbent[sel] # gov party
                        ,data$incumbent.emphasis.lag[sel]
                        ,data$prevvote[sel]
                        ,fes)
fakesweden2010 <- fakesweden2010[!is.na(coefs)]


prob_base <- prob_fake <- NULL
for(i in 1:nsim){
  # 1. predict true y (issue emphasis of conservative party 2010 in Sweden)
  # log odds ratios
  xb_base <- truesweden2010 %*% sims[i,]
  # get from log odds to probabilities
  odds_base <- exp(xb_base)
  prob_base[i] <- odds_base/(1+odds_base)
  
  # 2. predict fake y
  # log odds ratios
  xb_fake <- fakesweden2010 %*% sims[i,]
  # get from log odds to probabilities
  odds_fake <- exp(xb_fake)
  prob_fake[i] <- odds_fake/(1+odds_fake)
}

# define new dataset with fake issue emphasis
newdata <- data[cmp$date>201009 & cmp$parfam==60,]
newcmp <- cmp[cmp$date>201009 & cmp$parfam==60,]
length(unique(newcmp$country)) # 14 countries but 19 entries
# take the mean for more than one conservative party in system
# BUT: exclude elections in a polity that had already an election after 2010
# (Spain 201512, Estonia 201503, Latvia 201410)
newsel <- cmp$date>=201009 & cmp$parfam==60 & !c(cmp$countryname=="Spain" & cmp$date==201512) & !c(cmp$countryname=="Estonia" & cmp$date==201503) & !c(cmp$countryname=="Latvia" & cmp$date==201410) 
newdata <- data[newsel,]
newcmp <- cmp[newsel,]

newfes <- matrix(0,nrow=nrow(newdata),ncol=sum(grepl("^factor\\(cmp\\$",names(coefs))))
for(i in 1:nrow(newdata)){
  newfes[i,names(coefs[grepl("^factor\\(cmp\\$",names(coefs))])%in%paste0("factor(cmp$party)",newcmp$party[i])] <- 1
  newfes[i,names(coefs[grepl("^factor\\(cmp\\$",names(coefs))])%in%paste0("factor(cmp$year)",newcmp$year[i])] <- 1
}

newW.f.win <- W.f.win[newsel ,newsel]
newW.f.lose <- W.f.lose[newsel,newsel]
newW.d.win <- W.d.win[newsel,newsel]
newW.d.lose <- W.d.lose[newsel,newsel]

# Conservative parties with true input
trueconscenario <- cbind(1,newdata$issue.lag
                            ,newW.f.win %*% newdata$issue
                            ,newW.f.lose %*% newdata$issue
                            ,newW.d.win %*% newdata$issue
                            ,newW.d.lose %*% newdata$issue
                            ,newdata$vsgreen.lag
                            ,newdata$vsgreen.lag2
                            ,newdata$globalgreen.lag
                            ,newdata$gdp.lag
                            ,newdata$incumbent
                            ,newdata$incumbent.emphasis.lag
                            ,newdata$prevvote
                            ,newfes)
trueconscenario <- trueconscenario[,!is.na(coefs)]

# change actual issue emphasis for conservatives in Sweden to fake value
newdata$issue[newcmp$countryname=="Sweden" & newcmp$parfam==60 & newcmp$date==201009] <- median(prob_base)
# use these inputs to predict other countries (conservative parties)
fakeconscenario <- cbind(1,newdata$issue.lag
                            ,newW.f.win %*% newdata$issue
                            ,newW.f.lose %*% newdata$issue
                            ,newW.d.win %*% newdata$issue
                            ,newW.d.lose %*% newdata$issue
                            ,newdata$vsgreen.lag
                            ,newdata$vsgreen.lag2
                            ,newdata$globalgreen.lag
                            ,newdata$gdp.lag
                            ,newdata$incumbent
                            ,newdata$incumbent.emphasis.lag
                            ,newdata$prevvote
                            ,newfes)
fakeconscenario <- fakeconscenario[,!is.na(coefs)]

# 3. y with true and fake input
xb_all_true <- trueconscenario%*%coefs[!is.na(coefs)]
xb_all_fake <- fakeconscenario%*%coefs[!is.na(coefs)]
# get from log odds to probabilities (true and fake)
odds_all_true <- exp(xb_all_true)
odds_all_fake <- exp(xb_all_fake)
prob_all_true <- odds_all_true/(1+odds_all_true)
prob_all_fake <- odds_all_fake/(1+odds_all_fake)

tmp <- cbind(newcmp$country,newcmp$date,newcmp$party,prob_all_fake-prob_all_true)
colnames(tmp) <- c("country","date","party","change")
# remove sweden
tmp <- tmp[tmp[,"country"]!=11,]

# mean for each country
ch_country <- cbind(unique(tmp[,"country"]),NA)
for(i in unique(tmp[,"country"])){
  ch_country[match(i,unique(tmp[,"country"])),2] <- mean(tmp[tmp[,"country"]==i,"change"])
}
ch_country <- as.data.frame(ch_country
                            ,row.names=unique(newcmp$countryname[newcmp$country %in% ch_country[,1]]))
ch_country$region <- rownames(ch_country)
ch_country$region[ch_country$region=="United Kingdom"] <- "UK"


###################
# Figure 2
###################
world <- map_data("world")

tojoin <- as.data.frame(matrix(
  nrow = length(table(world$region)),
  ncol = 1,
  NA,
  dimnames = list(names(table(world$region)),NA)
))
tojoin$region <- rownames(tojoin)

all <- merge(tojoin,ch_country,by="region",all.x=T)
all <- all[,!is.na(colnames(all))]
colnames(all) <- c("region","country","issue")
all$issue <- all$issue*100

maps <- inner_join(world,all,by="region")
maps2 <- maps[maps$region %in% ch_country$region,]
countries <- unique(cmp$countryname)
countries[countries=="United Kingdom"] <- "UK"
maps3 <- maps[maps$region %in% countries,]

# PLOT
pdf("./Output/Figure_2.pdf",width=5,height=5)
draw <- ggplot(maps, aes(x = long, y = lat)) + theme(
  panel.background = element_rect(fill = "white",
                                  color = NA),
  panel.grid = element_blank(),
  axis.text.x = element_blank(),
  axis.text.y = element_blank(),
  axis.ticks = element_blank(),
  axis.title.x = element_blank(),
  axis.title.y = element_blank()
)

draw <- draw + coord_fixed(xlim = c(-9, 32.7),
                           ylim = c(36, 70),
                           ratio = 1.5)

draw <- draw + geom_polygon(data = maps,
                            aes(x = long,
                                y = lat,
                                group = group),
                            fill="white",
                            color = "black")

draw <- draw + geom_polygon(data = maps3,
                            aes(x = long,
                                y = lat,
                                group = group),
                            fill="grey",
                            color = "black")

draw <- draw + geom_polygon(data = maps2,
                            aes(fill = issue,
                                x = long,
                                y = lat,
                                group = group)
                            )
draw <- draw + theme(legend.position = "bottom"
                     ,legend.title = element_blank()
                     ,legend.text = element_text(size=8)
                     ,legend.key.width = unit(1, 'cm')
                     ,legend.key.height = unit(.3, 'cm')
                     ,legend.spacing=unit(.5, 'cm')
                     )
draw
dev.off()



###################
# Table B.1
###################
# alternative connectivity matrices based on ideological blocs
load("./Data/Matrices_bloc.RData")

# spatial lags
W.f.win <- W.f.win[complete,complete]
W.f.lose <- W.f.lose[complete,complete]
W.d.win <- W.d.win[complete,complete]
W.d.lose <- W.d.lose[complete,complete]
data$f.bloc.win <- W.f.win %*% data$issue
data$f.bloc.lose <- W.f.lose %*% data$issue
data$d.bloc.win.lag <- W.d.win %*% data$issue
data$d.bloc.lose.lag <- W.d.lose %*% data$issue

# Model 4
m4 <- glm(issue ~ issue.lag
          + f.bloc.win + f.bloc.lose
          + d.bloc.win.lag + d.bloc.lose.lag
          + vsgreen.lag
          + vsgreen.lag2
          + globalgreen.lag
          + gdp.lag
          + incumbent
          + incumbent.emphasis.lag
          + prevvote
          + factor(cmp$party)
          + factor(cmp$year)
          ,data=data
          ,family=quasibinomial("logit"))
coef4 <- sprintf("%.3f",round(coef(m4),3)[1:13])
se4 <- sprintf("%.3f",round(summary(m4)$coefficients[,"Std. Error"],3)[1:13])


# phi
phi <- c(round(summary(m3)$dispersion,4),round(summary(m4)$dispersion,4))
# Observations
obs <- c(length(m3$residuals),length(m4$residuals))
# RMSE
RMSE <- c(round(sqrt(mean((m3$fitted.values - m3$y)^2)),4)
          ,round(sqrt(mean((m4$fitted.values - m4$y)^2)),4))


cat("TABLE B.1: Table B.1: Fractional Logit Model Estimates of Environmental Issues Emphasis Based on Party Bloc Membership\n\n"
    ,append=F, file="./Output/TabB.1.txt")
cat("\t\t\t\t Model 3 \t\t Model 5\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("___________________________________________\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("Constant\t\t",coef3[1],"\t\t",coef4[1],"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("\t\t\t\t",paste0("(",se3[1],")"),"\t\t",paste0("(",se4[1],")"),"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("DV_t-1\t\t\t",coef3[2],"\t\t",coef4[2],"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("\t\t\t\t",paste0("(",se3[2],")"),"\t\t",paste0("(",se4[2],")"),"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("theta_1W^Fwy\t",coef3[3],"\t\t\t",coef4[3],"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("\t\t\t\t",paste0("(",se3[3],")"),"\t\t",paste0("(",se4[3],")"),"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("theta_2W^Fly\t",coef3[4],"\t\t",coef4[4],"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("\t\t\t\t",paste0("(",se3[4],")"),"\t\t",paste0("(",se4[4],")"),"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("theta_3W^Dwy\t",coef3[5],"\t\t\t",coef4[5],"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("\t\t\t\t",paste0("(",se3[5],")"),"\t\t",paste0("(",se4[5],")"),"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("theta_4W^Dly\t",coef3[6],"\t\t",coef4[6],"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("\t\t\t\t",paste0("(",se3[6],")"),"\t\t",paste0("(",se4[6],")"),"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("VSGreen_t-1\t\t",coef3[7],"\t\t",coef4[7],"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("\t\t\t\t",paste0("(",se3[7],")"),"\t\t",paste0("(",se4[7],")"),"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("VSGreen^2_t-1\t",coef3[8],"\t\t\t",coef4[8],"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("\t\t\t\t",paste0("(",se3[8],")"),"\t\t",paste0("(",se4[8],")"),"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("AvgGreen_n-1\t",coef3[9],"\t\t\t",coef4[9],"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("\t\t\t\t",paste0("(",se3[9],")"),"\t\t",paste0("(",se4[9],")"),"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("GDP_n-1\t\t\t",coef3[10],"\t\t\t",coef4[10],"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("\t\t\t\t",paste0("(",se3[10],")"),"\t\t",paste0("(",se4[10],")"),"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("Incumbent_t\t\t",coef3[11],"\t\t",coef4[11],"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("\t\t\t\t",paste0("(",se3[11],")"),"\t\t",paste0("(",se4[11],")"),"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("Incumbent_DV_t-1",coef3[12],"\t\t\t",coef4[12],"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("\t\t\t\t",paste0("(",se3[12],")"),"\t\t",paste0("(",se4[12],")"),"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("VS_t-1\t\t\t",coef3[13],"\t\t\t",coef4[13],"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("\t\t\t\t",paste0("(",se3[13],")"),"\t\t",paste0("(",se4[13],")"),"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("Party and\t\t","Yes \t\t\t Yes\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("Year FEs\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("___________________________________________\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("phi\t\t\t\t",phi[1],"\t\t",phi[2],"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("Observations\t",obs[1],"\t\t\t",obs[2],"\n"
    ,append=T, file="./Output/TabB.1.txt")
cat("RMSE\t\t\t",RMSE[1],"\t\t",RMSE[2]
    ,append=T, file="./Output/TabB.1.txt")


###################
# Table B.2
###################
### Average Marginal Effects
set.seed(123)
# Extract information for Simulations
coefs <- m3$coefficients[1:13]
vcov <- vcov(m3)[1:13,1:13]
nsim <- 10000
sims <- mvrnorm(nsim,mu=coefs,Sigma=vcov)
ame.sim <- matrix(NA,nrow=nsim,ncol=4)
colnames(ame.sim) <- c("f.bloc.win","f.bloc.lose"
                       ,"d.bloc.win.lag","d.bloc.lose.lag")
for(i in 1:nsim){
  # AME of spatial coefficients
  scenario <- cbind(1,data$issue.lag
                    ,data$f.bloc.win
                    ,data$f.bloc.lose
                    ,data$d.bloc.win.lag
                    ,data$d.bloc.lose.lag
                    ,data$vsgreen.lag,data$vsgreen.lag2
                    ,data$globalgreen.lag
                    ,data$gdp.lag,data$incumbent
                    ,data$incumbent.emphasis.lag
                    ,data$prevvote)
  # xb
  xb <- scenario %*% sims[i,]
  # get probabilities
  p <- exp(xb)/(1+exp(xb))
  
  # AMEs spatial short-term effects
  ame.sim[i,"f.bloc.win"] <- 1/nrow(data) * sum((sims[i,"f.bloc.win"]*s.f.bloc.win)*p*(1-p))
  ame.sim[i,"f.bloc.lose"] <- 1/nrow(data) * sum((sims[i,"f.bloc.lose"]*s.f.bloc.lose)*p*(1-p))
  ame.sim[i,"d.bloc.win.lag"] <- 1/nrow(data) * sum((sims[i,"d.bloc.win.lag"]*s.d.bloc.win.lag)*p*(1-p))
  ame.sim[i,"d.bloc.lose.lag"] <- 1/nrow(data) * sum((sims[i,"d.bloc.lose.lag"]*s.d.bloc.lose.lag)*p*(1-p))
}

ame <- matrix(NA,nrow=4,ncol=3)
colnames(ame) <- c("lowerCI","Estimate","upperCI")
rownames(ame) <- c("f.bloc.win","f.bloc.lose","d.bloc.win.lag","d.bloc.lose.lag")
ame[,"Estimate"] <- apply(ame.sim,2,mean)
ame[,"lowerCI"] <- apply(ame.sim,2,quantile, probs=.025)
ame[,"upperCI"] <- apply(ame.sim,2,quantile, probs=.975)

nams <- c("theta_1W^Fwy","theta_2W^Fly","theta_3W^Dwy","theta_4W^Dly")


cat("TABLE B.2: Average Marginal Effect Estimates of the Spatial Short-Term Effects (Ideo- logical Blocs)\n\n"
    ,append=F, file="./Output/TabB.2.txt")
cat("_____________________________________________\n"
    ,append=T, file="./Output/TabB.2.txt")
cat("\t\t\t\t Estimate \t\t 95%-CI\n"
    ,append=T, file="./Output/TabB.2.txt")
cat("_____________________________________________\n"
    ,append=T, file="./Output/TabB.2.txt")
cat(nams[1],"\t",sprintf("%.3f",ame[1,"Estimate"]),"\t\t"
    ,paste0("[",sprintf("%.3f",ame[1,"lowerCI"]),";\t",sprintf("%.3f",ame[1,"upperCI"]),"]")
    ,"\n"
    ,append=T, file="./Output/TabB.2.txt")
cat(nams[2],"\t",sprintf("%.3f",ame[2,"Estimate"]),"\t"
    ,paste0("[",sprintf("%.3f",ame[2,"lowerCI"]),";\t",sprintf("%.3f",ame[2,"upperCI"]),"]")
    ,"\n"
    ,append=T, file="./Output/TabB.2.txt")
cat(nams[3],"\t",sprintf("%.3f",ame[3,"Estimate"]),"\t\t"
    ,paste0("[",sprintf("%.3f",ame[3,"lowerCI"]),";\t",sprintf("%.3f",ame[3,"upperCI"]),"]")
    ,"\n"
    ,append=T, file="./Output/TabB.2.txt")
cat(nams[4],"\t",sprintf("%.3f",ame[4,"Estimate"]),"\t"
    ,paste0("[",sprintf("%.3f",ame[4,"lowerCI"]),";\t",sprintf("%.3f",ame[4,"upperCI"]),"]")
    ,"\n"
    ,append=T, file="./Output/TabB.2.txt")
cat("_____________________________________________\n"
    ,append=T, file="./Output/TabB.2.txt")


